function cholesky_decomposition()
% CHOLESKY_DECOMPOSITION Cholesky分解演示
% 
% 演示对称正定矩阵的Cholesky分解和LDL分解

fprintf('=== Cholesky分解和LDL分解 ===\n\n');

%% 1. Cholesky分解基础
fprintf('1. Cholesky分解基础\n');
fprintf('对于对称正定矩阵A，存在唯一的下三角矩阵L使得 A = L*L^T\n\n');

% 创建对称正定矩阵
A = [4, 2, 1; 
     2, 3, 0.5; 
     1, 0.5, 2];

fprintf('对称正定矩阵 A =\n');
disp(A);

% 检查对称性和正定性
fprintf('对称性检查: ||A - A^T|| = %e\n', norm(A - A', 'fro'));
eigenvals = eig(A);
fprintf('特征值: [%.4f, %.4f, %.4f]\n', eigenvals);
fprintf('最小特征值: %.4f (正定要求 > 0)\n', min(eigenvals));

% 执行Cholesky分解
L = cholesky_decomp(A);
fprintf('Cholesky因子 L =\n');
disp(L);

% 验证分解
LLT = L * L';
fprintf('验证 L*L^T - A 的误差: %e\n', norm(LLT - A, 'fro'));

%% 2. 与LU分解的比较
fprintf('\n2. 与LU分解的比较\n');
fprintf('对于对称正定矩阵，Cholesky分解比LU分解更高效\n\n');

% LU分解
[L_lu, U_lu] = lu_decomp_basic(A);
fprintf('LU分解:\n');
fprintf('L_LU =\n'); disp(L_lu);
fprintf('U_LU =\n'); disp(U_lu);

% 观察U = L^T的关系
fprintf('U_LU - L_LU^T 的误差: %e\n', norm(U_lu - L_lu', 'fro'));
fprintf('这验证了对称正定矩阵的LU分解中 U = L^T\n');

% 运算量比较
n = size(A, 1);
flops_cholesky = n^3/3;
flops_lu = 2*n^3/3;
fprintf('Cholesky分解运算量: %.0f\n', flops_cholesky);
fprintf('LU分解运算量: %.0f\n', flops_lu);
fprintf('Cholesky分解节省: %.1f%%\n', (1 - flops_cholesky/flops_lu) * 100);

%% 3. LDL分解
fprintf('\n3. LDL分解\n');
fprintf('对于对称矩阵A，存在 A = L*D*L^T\n');
fprintf('其中L是单位下三角矩阵，D是对角矩阵\n\n');

[L_ldl, D] = ldl_decomp(A);
fprintf('LDL分解:\n');
fprintf('L =\n'); disp(L_ldl);
fprintf('D =\n'); disp(D);

% 验证分解
LDLT = L_ldl * D * L_ldl';
fprintf('验证 L*D*L^T - A 的误差: %e\n', norm(LDLT - A, 'fro'));

% LDL与Cholesky的关系
D_sqrt = diag(sqrt(diag(D)));
L_from_ldl = L_ldl * D_sqrt;
fprintf('从LDL得到的Cholesky因子与直接计算的差异: %e\n', ...
        norm(L_from_ldl - L, 'fro'));

%% 4. 使用Cholesky分解求解线性方程组
fprintf('\n4. 使用Cholesky分解求解线性方程组\n');
fprintf('求解 A*x = b，分两步:\n');
fprintf('1) 求解 L*y = b (前向替换)\n');
fprintf('2) 求解 L^T*x = y (后向替换)\n\n');

b = [7; 5.5; 3.5];
fprintf('右端向量 b = [%g; %g; %g]\n', b);

% 第一步：前向替换 L*y = b
y = forward_substitution(L, b);
fprintf('中间解 y = [%.6f; %.6f; %.6f]\n', y);

% 第二步：后向替换 L^T*x = y
x = backward_substitution(L', y);
fprintf('最终解 x = [%.6f; %.6f; %.6f]\n', x);

% 验证解
residual = A * x - b;
fprintf('残差 ||A*x - b|| = %e\n', norm(residual));

% 与MATLAB内置求解器比较
x_matlab = A \ b;
fprintf('MATLAB解: x = [%.6f; %.6f; %.6f]\n', x_matlab);
fprintf('差异: ||x_chol - x_matlab|| = %e\n', norm(x - x_matlab));

%% 5. 修正Cholesky分解
fprintf('\n5. 修正Cholesky分解\n');
fprintf('当矩阵不是正定时，可以使用修正Cholesky分解\n\n');

% 创建不定矩阵
A_indef = [1, 2, 3; 
           2, 1, 4; 
           3, 4, 1];

fprintf('不定矩阵 A =\n');
disp(A_indef);
eigenvals_indef = eig(A_indef);
fprintf('特征值: [%.4f, %.4f, %.4f]\n', eigenvals_indef);

% 尝试标准Cholesky分解
try
    L_fail = cholesky_decomp(A_indef);
    fprintf('标准Cholesky分解成功\n');
catch ME
    fprintf('标准Cholesky分解失败: %s\n', ME.message);
end

% 修正Cholesky分解
[L_modified, E] = modified_cholesky(A_indef);
fprintf('修正Cholesky分解:\n');
fprintf('L =\n'); disp(L_modified);
fprintf('修正矩阵 E =\n'); disp(E);

% 验证 A + E = L*L^T
A_modified = A_indef + E;
LLT_modified = L_modified * L_modified';
fprintf('验证 (A+E) - L*L^T 的误差: %e\n', norm(A_modified - LLT_modified, 'fro'));
fprintf('修正矩阵的范数: %e\n', norm(E, 'fro'));

%% 6. 性能测试
fprintf('\n6. 性能测试\n');
fprintf('比较不同分解方法的性能\n\n');

sizes = [100, 200, 400, 800];
fprintf('矩阵规模    Cholesky    LDL分解    LU分解\n');
fprintf('--------    --------    -------    ------\n');

for i = 1:length(sizes)
    n = sizes(i);
    
    % 生成随机对称正定矩阵
    A_test = randn(n, n);
    A_test = A_test' * A_test + eye(n);  % 确保正定
    
    % Cholesky分解
    tic;
    L_test = cholesky_decomp(A_test);
    time_chol = toc;
    
    % LDL分解
    tic;
    [L_ldl_test, D_test] = ldl_decomp(A_test);
    time_ldl = toc;
    
    % LU分解
    tic;
    [L_lu_test, U_lu_test] = lu_decomp_basic(A_test);
    time_lu = toc;
    
    fprintf('%6d      %6.4f s    %6.4f s    %6.4f s\n', ...
            n, time_chol, time_ldl, time_lu);
end

%% 7. 数值稳定性
fprintf('\n7. 数值稳定性\n');
fprintf('Cholesky分解的数值稳定性分析\n\n');

% 创建病态对称正定矩阵
n = 5;
D_ill = diag(logspace(-10, 0, n));  % 条件数很大
Q = orth(randn(n, n));  % 正交矩阵
A_ill = Q * D_ill * Q';
A_ill = (A_ill + A_ill') / 2;  % 确保对称

fprintf('病态矩阵条件数: %.2e\n', cond(A_ill));

% Cholesky分解
L_ill = cholesky_decomp(A_ill);
LLT_ill = L_ill * L_ill';
chol_error = norm(LLT_ill - A_ill, 'fro') / norm(A_ill, 'fro');

fprintf('Cholesky分解相对误差: %e\n', chol_error);

% 求解线性方程组
b_ill = randn(n, 1);
y_ill = forward_substitution(L_ill, b_ill);
x_ill = backward_substitution(L_ill', y_ill);

residual_ill = A_ill * x_ill - b_ill;
rel_residual = norm(residual_ill) / (norm(A_ill, 'fro') * norm(x_ill));

fprintf('求解相对残差: %e\n', rel_residual);

end

function L = cholesky_decomp(A)
% CHOLESKY_DECOMP Cholesky分解实现
% 输入: A - 对称正定矩阵
% 输出: L - 下三角Cholesky因子

n = size(A, 1);
L = zeros(n, n);

for k = 1:n
    % 计算对角元素
    sum_sq = sum(L(k, 1:k-1).^2);
    if A(k, k) - sum_sq <= 0
        error('矩阵不是正定的');
    end
    L(k, k) = sqrt(A(k, k) - sum_sq);
    
    % 计算第k列的下三角部分
    for i = k+1:n
        sum_prod = sum(L(i, 1:k-1) .* L(k, 1:k-1));
        L(i, k) = (A(i, k) - sum_prod) / L(k, k);
    end
end

end

function [L, D] = ldl_decomp(A)
% LDL_DECOMP LDL分解实现
% 输入: A - 对称矩阵
% 输出: L - 单位下三角矩阵, D - 对角矩阵

n = size(A, 1);
L = eye(n);
D = zeros(n, n);

for k = 1:n
    % 计算对角元素
    sum_term = sum(L(k, 1:k-1).^2 .* diag(D(1:k-1, 1:k-1))');
    D(k, k) = A(k, k) - sum_term;
    
    % 计算第k列的下三角部分
    for i = k+1:n
        sum_term = sum(L(i, 1:k-1) .* L(k, 1:k-1) .* diag(D(1:k-1, 1:k-1))');
        L(i, k) = (A(i, k) - sum_term) / D(k, k);
    end
end

end

function [L, E] = modified_cholesky(A)
% MODIFIED_CHOLESKY 修正Cholesky分解
% 对于不定矩阵，找到最小的修正E使得A+E正定

n = size(A, 1);
L = zeros(n, n);
E = zeros(n, n);
gamma = 0;
xi = 0;

% 计算参数
for i = 1:n
    gamma = max(gamma, abs(A(i, i)));
    for j = i+1:n
        xi = max(xi, abs(A(i, j)));
    end
end

delta = max(eps, xi^2 / gamma);
beta = sqrt(max(gamma, xi / sqrt(n^2 - 1), eps));

for k = 1:n
    % 计算对角元素的修正
    sum_sq = sum(L(k, 1:k-1).^2);
    theta = 0;
    for i = k+1:n
        sum_prod = sum(L(i, 1:k-1) .* L(k, 1:k-1));
        c_ik = A(i, k) - sum_prod;
        theta = max(theta, abs(c_ik));
    end
    
    diagonal_val = A(k, k) - sum_sq;
    modification = max(0, -diagonal_val + max(theta^2 / beta, delta));
    
    if modification > 0
        E(k, k) = modification;
    end
    
    L(k, k) = sqrt(diagonal_val + modification);
    
    % 计算第k列的下三角部分
    for i = k+1:n
        sum_prod = sum(L(i, 1:k-1) .* L(k, 1:k-1));
        L(i, k) = (A(i, k) - sum_prod) / L(k, k);
    end
end

end

function [L, U] = lu_decomp_basic(A)
% 基本LU分解
n = size(A, 1);
L = eye(n);
U = A;

for k = 1:n-1
    for i = k+1:n
        L(i, k) = U(i, k) / U(k, k);
        U(i, k:n) = U(i, k:n) - L(i, k) * U(k, k:n);
    end
end
end

function x = forward_substitution(L, b)
% 前向替换
n = size(L, 1);
x = zeros(n, 1);
for i = 1:n
    x(i) = (b(i) - L(i, 1:i-1) * x(1:i-1)) / L(i, i);
end
end

function x = backward_substitution(U, b)
% 后向替换
n = size(U, 1);
x = zeros(n, 1);
for i = n:-1:1
    x(i) = (b(i) - U(i, i+1:n) * x(i+1:n)) / U(i, i);
end
end